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SECTION  1 
INTRODUCTION 


The  calculation  of  the  large  deformation  response  of  structures  by  finite  element 
methods  is  often  a  difficult  task  because  the  user  has  little  guidance  as  to  the  size  of  the 
elements  that  need  to  be  used.  This  is  a  particularly  severe  handicap  in  calculating  the 
failure  of  structures  because  failtue  often  results  in  localization  of  deformation. 
Localization  requires  considerable  refinement  of  the  mesh,  in  the  area  of  localization  if 
the  failure  mode  is  to  be  calculated  accurately.  Therefore,  the  user  must  anticipate 
where  localization  is  likely  to  occur  in  designing  the  mesh  or  make  several  runs  to 
achieve  an  accurate  solution.  Adaptivity  provides  a  way  for  reducing  this  burden  by 
providing  an  automatic  selection  of  appropriate  element  refinement.  In  adaptive 
programs,  elements  are  subdivided  in  areas  where  refinement  is  needed  so  that  the 
required  accuracy  is  maintained.  This  refinement  is  guided  by  error  indicators,  which 
indicate  an  approximate  level  of  the  errors  which  are  engendered  by  the  finite  element 
approximation;  refinement  criteria  which  are  not  based  on  error  criteria  can  also  be 
used. 

This  report  describes  the  implementation  of  adaptivity  in  the  computer  program 
DYNA2D,  (Whir ley,  et  al  1992).  The  adaptivity  is  performed  by  subdividing  elements 
into  smaller  elements  according  to  error  or  refinement  criteria,  which  is  called  h- 
adaptivity.  The  name  /j-adaptivity  stems  from  the  association  of  h  with  element  size  in 
the  finite  element  literature.  DYNA2D  is  a  program  originated  by  John  Hallquist  at 
Lawrence  Livermore  Laboratory.  This  program  is  a  general  purpose  program  for  the 
nonlinear  dynamics  of  solids  and  structures.  It  is  based  on  explicit  time  integration,  so 
it  does  not  store  or  compute  a  stiffness  matrix.  H-adaptivity  enables  the  user  of  the 
program  to  automatically  refine  portions  of  the  mesh  according  to  error  or  refinement 
criteria  so  that  accuracy  is  increased  for  a  given  amoimt  of  computational  resources. 

This  is  of  particular  importance  in  nonlinear  failure  problems  because  the  most  crucial 
areas  for  the  placement  of  nodes  and  refinement  of  the  mesh  are  in  the  failure  areas.  In 
li-adaptivity,  the  program  will  automatically  refine  the  mesh  in  areas  of  the  mesh  in 
which  the  estimated  error  is  large,  and  thus  improve  accuracy. 

The  program  DYNA2D  was  chosen  for  this  effort  because  it  is  widely  used  in  the 
defense  commvmity  and  has  a  large  range  of  capabilities  for  nonlinear  analysis  of 
weapons  effects  and  related  problems.  Among  the  attractive  features  and  capabilities  of 
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DYNA2D  are  a  large  material  library,  an  effective  sliding  interface  for  contact-impact, 
and  the  availability  of  capabilities  such  as  rezoning  and  remeshing,  postprocessing  and 
interactive  display,  and  the  availability  of  a  compatible  preprocessor  MAZE,  (Hallquist 
1983).  Furthermore,  the  development  of  these  procedures  in  DYNA2D  has  provided 
insight  into  adding  adaptivity  to  other  programs  such  as  DYNA2D. 

During  the  course  of  this  work,  it  was  foimd  that  it  is  very  difficult  to  modify  a  program 
such  as  DYNA2D  to  internally  implement  /i-adaptivity.  The  implementation  of 
adaptivity  requires  so  many  changes  throughout  the  program  that  our  initial  efforts  in 
this  direction  were  totally  unsuccessful  and  could  not  even  invoke  a  very  limited  subset 
of  DYNA2D  functionality.  After  further  consideration  of  this  approach,  it  also  became 
clear  that  internal  modifications  would  preclude  its  use  by  many  other  users,  since 
many  versions  of  DYNA2D  are  imder  separate  development  in  the  defense  community. 

Therefore,  we  selected  an  approach  in  which  the  programming  for  the  /i-adaptivity  is 
completely  external  to  DYNA2D  and  an  li-adaptive  nm  is  executed  through  a  series  of 
DYNA2D  restarts  and  rezones.  This  enabled  us  to  keep  all  coding  for  die  /i-adaptivity 
completely  separate  from  DYNA2D,  thus  making  the  /i-adaptive  procedures  usable  in 
later  versions  of  DYNA2D  and  versions  of  DYNA2D  which  have  been  modified  by 
other  agencies. 

This  report  is  organized  as  follows.  Section  2  describes  the  general  framework  of 
DYNA2D,  of  /i-adaptivity  and  the  procedures  involved  in  its  implementations  and  the 
error  criteria  which  have  been  implemented  in  the  program.  Section  3  gives  some 
examples  of  /i-adaptive  solutions,  which  is  followed  by  conclusions  in  Section  4. 
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SECTION  2 
METHODOLOGY 


2.1  OUTLINE  AND  NOTATION. 

This  Section  briefly  describes  the  methodology  of  DYNA2D  and  then  gives  a 
description  of  Ii-adaptivity,  the  reasons  for  selection  of  /z-adaptivity  and  the  refinement 
and  error  criteria  implemented  in  this  work. 

Vectors,  tensors  and  matrices  in  this  report  are  indicated  by  boldface.  Indicial  notation 
is  also  used.  In  indicial  notation,  the  elements  of  a  vector  tensor  or  matrix  are  indicated 
by  subscripts;  for  example,  in  indicial  notation  the  elements  of  vector  x  are  indicated  by 
X,  .  Upper  case  subscripts  are  used  for  nodal  numbers,  so  x-,  are  the  coordinates  at  node 
L 

2.2  COMPUTAHONAL  PROCEDURE  IN  DYNA2D. 

DYNA2D  is  a  finite  element  program  for  calculating  the  nonlinear  response  of  solids  in 
two  dimensions;  in  addition,  axisymmetric  three  dimensional  structures  can  be  treated. 
The  program  uses  explicit  time  integration.  In  the  nonlinear  continuum  mechanics 
formulation  in  DYNA2D,  the  deformation  of  the  solid  is  described  by 

X  =  x(X,  t)  (2.1) 

where  x=[x,y]  are  the  current  spatial  coordinates  and  X=[X,Y]  are  the  material  or 
original  coordinates  of  a  material  point.  In  a  finite  element  program,  the  deformation  is 
approximated  by  subdividing  the  domain  of  the  problem  Q  into  elements  such  that 
Q  =  In  each  element,  the  deformation  is  approximated  by  shape  functions  N,  so 
that 


/=! 

where  rj  are  the  coordinates  of  the  element  in  the  parent  element  plane  and  Xy  are  the 
coordinates  of  node  I  of  element  e:  x,,  =  [x/,  y‘  ];  m  is  the  number  of  nodes  in  the  element. 


The  displacement  field  u(X,t)  is  approximated  by  the  same  shape  functions 
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(2.3) 


/=! 

where  are  tiie  nodal  displacements  of  element  e;  u-,  =  The  nodal  velocities 

are  obtained  from  Equation  (2.3)  by  talcing  a  time  derivative,  which  gives 


v'(X.  t)  =  X  77)v' (0  =  (2.4) 

7=1  7=1 

where  are  the  nodal  velocities  of  element  e.  The  nodal  displacements  and  nodal 
velocities  of  the  entire  mesh  are  denoted  by  Ug  and  v^,  respectively;  in  vector  notation, 
they  are  denoted  by  and  y,,  respectively. 

The  deformation  in  DYNA2D  is  measured  by  the  velocity  strain  tensor,  which  is 
denoted  by  e  and  given  by 


1 


dxj  dXf 


or  £  =  synNy 


(2.5) 


The  velocity  strain  is  a  rate  of  measure  of  deformation,  rather  than  a  measure  of  the  total 
deformation,  or  strain,  which  has  occurred  at  a  point.  The  rate-of-deformation  in  an 
element  given  by 


m  rn 

e  =  =  ^B,y, 


7=1 


7=1 


(2.6) 


where 


B;  =  symVN,  (2.7) 

In  explicit  programs,  rate  measures  of  deformation  at  times  oscillate  markedly  in 
comparison  to  total  measures  of  deformation.  Furthermore,  the  integral  of  the  velocity 
strain  is  not  path  independent,  so  the  integral  of  the  rate-of-deformation  does  not 
represent  a  good  measure  of  the  total  strain. 

The  finite  element  procedure  in  DYNA2D  uses  a  semidiscretization  of  the  momentum 
equation.  This  leads  to  the  finite  element  form  of  the  equations  of  motion 

M,v,=fr-f;  =f,  (2.8) 
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where  M,  is  a  mass  matrix,  which  is  lumped,  and  consequently  diagonal;  a  superposed 
dot  denotes  a  time  derivative  and  f  and  f  are  the  external  and  internal  nodal  forces; 
as  indicated  in  the  second  equality,  the  difference  between  the  external  and  internal 
nodal  forces  are  often  indicated  as  a  net  force.  The  internal  nodal  forces  are  given  by 

tf  =  {2.9) 

o. 

where  a  is  the  Cauchy  stress. 

Equation  (2.8)  is  integrated  in  time  by  a  central-difference  explicit  method.  The  central 
difference  metiiod  is  conditionally  stable,  and  its  time  step  depends  on  the  smallest 
element  in  the  mesh.  Most  simulations  of  nonlinear  structural  dynamics  problems 
require  many  time  steps.  Generally,  on  the  order  of  10^  or  10^  time  steps  are  used,  but 
it  is  not  imusual  to  have  simulations  with  10*  time  steps.  This  has  important 
implications  on  the  implementation  of  adaptivity,  since  the  mesh  can  be  refined  only  on 
the  order  of  ten  times  in  a  simulation  if  the  computational  cost  of  adaptivity  is  to  be 
reasonable. 

2.3  ADAPTIVITY. 

We  will  first  describe  why  /i-adaptivity  was  selected  for  DYNA2D  among  the  various 
types  of  adaptivity  available  for  finite  elements  analysis.  There  are  basically  three  types 
of  adaptivity: 

1.  r-adaptivity,  in  which  the  nodes  are  moved  so  that  the  degrees  of  freedom  are 
concentrated  in  the  areas  where  they  are  needed; 

2.  p-adaptivity,  in  which  the  order  of  the  interpolating  function  in  the  element  is 
increased  when  more  accuracy  is  needed; 

3.  Iz-adaptivity,  in  which  the  elements  are  subdivided  when  more  accuracy  is  needed. 

The  capabilities  of  r-adaptivity  are  quite  limited  in  most  practical  meshes.  For  example, 
if  the  nodes  are  needed  primarily  on  a  diagonal  line  across  a  square  mesh,  when  the 
starting  point  is  a  uniform  mesh,  it  is  impossible  to  move  many  of  the  nodes  to  the 
diagonal  without  inducing  excessive  mesh  distortion.  Mesh  distortion  degrades  the 
accuracy  of  finite  elements,  so  some  of  the  benefits  of  moving  the  nodes  are  lost. 
Furthermore,  nodes  cannot  be  moved  across  material  interfaces  without  invalidating 
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the  model.  Therefore,  in  general,  very  little  additional  accuracy  can  be  gained  by  r- 
adaptivity. 

In  p-adaptivity,  the  order  of  the  interpolating  polynomial  is  raised  where  more  accuracy 
is  needed.  For  example  if  an  element  is  initially  a  bilinear  element,  then  p-adaptivity 
would  increase  the  order  of  the  element  to  a  biquadratic  element.  P-adaptivity  is 
unsuitable  for  this  effort  with  explicit  nonlinear-structural  dynamics  for  two  reasons: 

1.  It  is  not  compatible  with  for  the  approach  chosen  for  the  implementation  of 
adaptivity,  namely  a  paradigm  in  which  all  modifications  for  adaptivity  are  completely 
outside  of  the  target  computer  program,  because  DYNA2D  contains  only  low  order 
elements. 

2.  This  type  of  adaptivity  is  not  advisable  in  explicit  codes  because  explicit  methods  use 
diagonal  mass  matrices,  and  diagonal  mass  matrices  of  suitable  accuracy  have  not  been 
developed  for  higher  order  elements. 

Based  on  these  arguments,  /j-adaptivity  was  selected  for  this  implementation  of 
DYNA2D.  In  li-adaptivity,  it  is  also  possible  to  fuse  elements,  or  in  other  words,  to 
xmrefine  the  mesh,  when  a  coarser  mesh  is  adequate.  However,  only  refinement,  which 
we  will  often  call  fission,  is  included  in  this  implementation  of  /i-adaptivity.  The 
elements  are  not  fused  when  the  refinement  is  no  longer  necessary,  because  fusion  of 
elements  requires  a  complex  data  base  and  is  quite  difficult  for  general  purpose 
programs. 

The  fission  process  for  the  elements  in  DYNA2D  is  shown  in  Figure  2-1.  As  can  be  seen 
from  the  figure,  triangles  are  refined  by  splitting  them  into  four  triangles  by  connecting 
the  midsides  of  the  nodes.  Quadrilaterals  are  refined  by  splitting  them  into  four 
quadrilateral  by  connecting  the  midpoints  of  the  sides.  In  the  fission  process,  new 
nodes  are  also  generated.  As  can  be  seen  from  Figure  2-1,  the  nodes  generated  by 
refinement  are  either  free  nodes  or  slave  nodes.  The  motion  of  free  nodes  is  completely 
xmconstrained  except  for  boundary  conditions,  so  the  addition  of  these  nodes  adds  extra 
degrees  of  freedom  to  the  model.  The  motion  of  slave  nodes  is  constrained  by  the  need 
to  observe  compatibility  between  adjacent  elements.  The  nodes  at  the  midpoints  of 
sides  of  the  larger  elements  are  slave  nodes;  the  treatment  of  these  nodes  is  described 
later. 
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Each  subdivision  is  called  a  level  of  fission.  Elements  may  be  refined  an  arbitrary 
number  of  levels,  so  that  elements  which  are  formed  by  splitting  or  subdivision  of  a 
previous  element  may  be  further  subdivided.  The  number  of  subdivisions  is  limited  in 
the  program  by  an  input  parameter. 

The  fission  process,  the  arrangement  of  elements  is  limited  by  the  1-irregular  rule 
(Devloo,  et  al  1987).  This  rule  requires  that  the  difference  in  the  level  of  subdivision 
between  adjacent  elements  be  at  most  two.  As  a  consequence,  at  most  one  slave  node 
can  be  inserted  in  any  edge. 


The  motion  of  slave  nodes  is  constrained  by  the  motion  of  the  adjacent  comer  nodes  so 
that  compatibility  is  maintained  between  the  subdivided  element  and  the  larger  element 
adjacent  to  it.  The  displacement  and  velocity  fields  for  elements  in  DYNA2D,  the  three 
node  triangle  and  four  node  quadrilateral,  are  linear  along  the  edges.  The  motion  of  a 
slave  node  must  conform  to  the  linear  velocity  field  of  the  larger  element  to  satisfy 
compatibility.  For  example,  if  we  consider  nodes  1, 2  and  3  in  Figure  2-1,  then  to 
maintain  compatibility  of  the  adjacent  elements  along  this  edge  the  velocity  of  the 
midpoint  node,  that  is  the  slave  node,  must  be  given  by 

'^2=^(Vi+V3)  (2.10) 

where  Vj  are  the  nodal  velocities;  =  [v^/,  v,.,  j  where  u*;  and  v^i  are  the  x  and  y 

components  of  nodes  I,  respectively.  The  equations  of  motion.  Equations  2.8,  do  not 
pertain  to  slave  nodes.  Therefore,  the  forces  on  the  slave  nodes  must  be  transferred  to 
master  nodes.  Nodes  which  are  not  constrained,  such  as  nodes  1  and  3,  are  called 
master  nodes.  The  internal  force  at  node  2  is  transferred  to  nodes  1  and  3  according  to 
the  following  procedure  bcised  on  the  principle  of  virtual  work.  Let  the  internal  forces 
at  node  I  be  denoted  by  f  f  f '^'j.  By  the  principle  of  virtual  work,  the  virtual  work 

of  nodal  force  transferred  from  slave  node  2  to  the  virtual  work  of  the  nodal  forces 
transferred  to  master  nodes  1  and  3.  This  means  that 

Sy^  ■  f f  =  Sv,  •  f f'  +  5v3  •  f f  (2.11) 
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From  Equation  (2.10), 

S\2  =  +  5V3)  (2.12) 

so  substituting  the  above  into  Equation  (2.11)  yields 

-(5v,  +  5V3)  •  ff  =  5v,  •  ff'  +  5\,  •  (2.13) 

2 

The  above  must  hold  for  arbitrary  5v,  and  5V3  so  it  follows  that 

f,«  ^  1  (2.14) 

'  2  ^ 

(2.15) 

2 

These  nodal  forces  are  added  to  the  nodal  forces  which  arise  at  nodes  1  and  from  the 
elements  directly  cormected  to  these  nodes. 

®  slave  node  O  free  node 


Figure  2-1.  Example  of  ^-adaptivity  for  a  quadrilateral  element;  the  left  element  has 
been  refined  one  level. 
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2.4  ERROR  AND  REFINEMENTINDICATORS. 


In  adaptive  methods,  an  error  or  refinement  indicator  must  be  computed  to  determine 
which  elements  need  to  be  refined.  Error  indicators  are  distinguished  from  refinement 
indicators  in  that  the  latter  make  no  attempt  to  estimate  an  error  but  simply  refine  on 
the  basis  of  certain  characteristics  of  the  solution.  In  nonlinear,  transient  analysis,  little 
theoretical  justification  is  available  for  the  error  indicators  which  are  used.  However, 
preliminary  studies  (Belytschko  and  Tabbara  1993)  show  that  the  indicators  used  in 
linear  analysis  are  also  quite  effective  in  nonlinear  analysis.  In  addition,  to  an  error 
indicator,  several  refinement  indicators  are  included  in  this  program.  These  refinement 
indicators  use  the  magnitude  of  specific  physical  quantities  to  drive  the  adaptivity. 

The  error  indicator  selected  here  is  based  on  an  estimate  of  the  error  in  the  Green- 
Lagrange  strain.  The  estimate  is  obtained  by  a  projection  of  the  finite  element  solution 
with  the  moving  least-square  interpolation  described  in  (Tabbara,  et  al  1994).  The 
projection  criterion  developed  here  is  very  similar  to  the  local  projection  error  indicator 
(Zienkiewicz  and  Zhu  1987)  except  that  the  projected  solution  is  obtained  by  a  moving 
least-square  approximation,  which  has  been  found  to  be  a  more  accurate  in  (Tabbara,  et 
al  1994).  This  error  criterion  is  an  improvement  of  the  earlier  global  projections  of 
(Zienkiewicz  and  Zhu  1987).  This  procedure  is  also  known  as  error  estimation  by  a 
recovery  technique. 

In  a  recovery  technique,  a  local  projection  is  used  to  obtain  an  approximation  to  the 
exact  solution  denoted  by  u*(X,t).  This  projected  solution  is  obtained  by  minimizing  the 
L2  norm  of  the  difference  between  the  projected  i.e.  recovered  solution  and  the  finite 
element  solution.  The  error  in  the  Green-Lagrange  stain  in  an  element  e  of  any  function 
f(X,t)  is  denoted  by  6^,  and  is  given  by 

(2.16) 

a. 

where  E*j  is  obtained  from  the  projected  solution  and  is  the  Green-strain  computed 
from  the  finite  element  solution.  As  can  be  seen  from  the  above,  die  integral  of  the  error 
is  divided  by  the  domain  (or  area  in  two  dimensions)  of  the  element,  so  it  is  a  density  of 
error. 
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The  rate-of-deformation  tensor  which  is  used  in  DYNA2D  for  constitutive  evaluations 
Equation  (2.5)  is  not  suitable  as  an  error  criterion  because  it  is  not  path-independent, 
furthermore,  since  it  gives  a  rate  of  strain  it  tends  to  be  noisy.  The  Green-Lagrange 
strain  is  path  independent  and  a  measure  of  strain.  Therefore,  it  is  more  suitable  for  our 
own  error  criterion.  The  Green-Lagrange  strain  tensor  E  is  given  by 


du,  dUj  dut  du^ 

dx,  dx,  dx, 


dXj 


(2.17) 


where  is  the  deformation  gradient,  which  is  given  by 


dXj  _  dUj  g 

dx.  ■  dXj  " 


(2.18) 


and  6ij  is  the  Kronecker  delta,  or  unit  matrix. 

The  projected  Green-Lagrange  strain  is.  obtained  from  a  moving  least-square 
approximation  for  the  displacement  field.  For  the  pvupose  of  computing  this 
approximation,  a  local  circular  domain  Qg  is  defined  about  the  centroid  of  the  element 
of  interest.  The  radius  of  the  domain  is  denoted  by  it  is  set  by  a  procedure  described 
subsequently.  This  domain  is  called  the  domain  of  influence  of  the  approximation. 


The  projection  of  the  displacement  field  u*(X,  t)  is  obtained  by  a  weight  least  square 
projection.  For  this  purpose,  the  displacement  fields  are  approximated  in  the 
subdomain  £2^  by 

M;(X,0  =  P(X)a,(X,0  (2.19) 

where  P(X)  is  a  monomial  basis  and  a,.(X,t)  are  coefficients  which  vary  in  space.  A 
bilinear  basis  is  used  for  P(X) ,  so  P(X)  is  given  by 

P(X)  =  [l,XK.Xy]  (2.20) 

In  the  moving  least-square  approximation,  the  parameters  a,  (X,  t)  at  any  time  t  are 
determined  by  minimizing  the  weighted  quadratic  form 

5{a,(X,l))=  X“'(X-X,)(P(X,)a,(X,r)-«;,f  (2.21) 

/=1 
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where  are  the  components  of  the  finite  element  nodal  displacements  at  node  I  and 
^(X  -  X/)  is  a  weight  function  given  by 


»>{X-  X,)  =  e^HIX  -  X,|/(0.4<i.)’) 


(2.22) 


and  is  twice  the  maximum  element  diagonal  for  any  element  in  the  domain  of 
influence.  This  choice  for  is  based  on  experience  which  shows  that  the  domain  of 
influence  should  include  the  elements  which  surroimd  the  element  under  consideration 
but  that  the  exponential  should  have  decayed  markedly  at  the  nodes  which  are  not  part 
of  the  element  under  consideration.  The  minimization  of  5(a,  (X,  r))  leads  to  a  system  of 
linear  equations  in  a,  (X,  t).  The  gradient  of  the  displacement  field  is  then  obtained 
using  Equations  (2.19)  and  (2.20)  which  gives 


du 


dX 


(2.23) 


^34  +  ^4*-^ 


a2,,  +  a4,T 


(2.24) 

(2.25) 


du* 

-^  =  a3v  +  a4;,y 


(2.26) 


These  gradients  are  then  substituted  into  Equation  (2.17)  to  obtain  the  projected  Green- 
strain  E*y  Note  that  the  derivatives  of  the  coefficients  a,  are  not  considered  in 

computing  the  gradients. 


The  remaining  indicators  included  in  the  program  are  not  based  on  an  estimate  of  error 
but  are  refinement  indicators  which  drive  refinement  on  the  basis  of  certain  physical 
quantities  that  are  likely  to  be  high  where  refinement  is  needed.  The  following  are 
included: 


1.  The  effective  stress  which  is  given  by 


where  s-^  is  the  deviatoric  stress  which  is  given  by 


(2.27) 


11 


Sij  =  a,j-p6y 


(2.28) 


2.  The  trace  of  the  stress  tensor  which  is  given  by 

trace(ay)  =  (J,^  (2.29) 

3.  The  effective  plastic  strain  is  given  by 

e^=\’de'’  (2.30) 

Jo 

where 


(2 
-de^def^ 

y  w 


U 


(2.31) 


2.5  IMPLEMENTATION. 

H-adaptivity  was  implemented  in  DYNA2D  by  developing  an  external  program  which 
performs  the  calculation  of  error  criteria  and  the  remeshing.  The  details  of  the 
implementation  are  described  subsequently,  after  the  procedure  is  described. 

The  basic  procedure  of  ft-adaptivity  for  nonlinear  transient  problems  is  summarized  in 
Table  2-1.  As  can  be  seen  from  Table  2-1,  the  adaptive  program  updates  the  mesh  by 
integrating  the  equations  of  motion  a  selected  time  interval,  denoted  by  and  then 

checks  the  refinement  or  error  criteria  for  all  of  the  elements  in  the  mesh.  Any  elements 
for  which  the  refinement  criteria  are  exceeded  are  then  subdivided. 

The  program  then  repeats  the  time  interval  with  the  refined  mesh  and  writes  a 
restart  tape;  it  then  continues  for  the  time  interval  beyond  the  time  at  which  the 
error  criterion  was  previously  checked  before  checking  refinement  criteria  again.  If  the 
refinement  criteria  indicate  that  any  elements  need  to  be  refined,  the  elements  are 
subdivided  and  the  computation  is  restarted  with  the  restart  tape  which  was  previously 
written.  The  procedure  is  called  a  go-back  procedure. 
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Table  2-1.  Flowchart  of  Adaptive  Procedure. 

1.  Initialize,  set  t^,  = 

2.  Integrate  the  equations  of  motion.  Equations  2-7,  to  a  specified  time  . 

3.  Qieck  error  or  refinement  indicators  0,  for  all  elements  e, 

if  0^>refinement  threshold  for  element  e,  subdivide  element  e. 

4.  Set  +  At^,  and  go  to  2. 

A  schematic  of  the  go-back  procedure  is  shown  in  Figure  2-2,  which  shows  the  time  line 
for  a  computation  with  two  refinement  steps  involving  a  total  of  three  meshes:  Mesh  1, 
Mesh  2,  and  Mesh  3.  As  can  be  seen  from  Figme  2-2,  the  computation  starts  with  Mesh 
1,  which  is  refined  at  the  first  terror  check  to  obtain  Mesh  2;  Mesh  2  repeats  the 
computation  over  the  time  interval  from  (o  <  t  <  At„^^,^  and  then  continues  to  2At„^„p,. 
Time  histories  are  stored  for  the  interval  0  <  r  <  At^,  only  for  the  Mesh  2  computations, 
not  for  the  Mesh  1  computation.  At  2At^,,  the  next  error  check  is  made.  In  the 
literature,  for  example  (Devloo,  Oden  and  Strouboulis  1987),  it  is  often  advocated  that 
the  computation  not  proceed  beyond  the  last  error  check  tmtil  the  error  criteria  are 
satisfied.  Unfortimately,  this  leads  to  the  possibility  that  the  computation  will  iterate 
many  times  in  a  single  time  interval,  which  is  not  desirable  in  a  production  program. 


Figure  2-2.  Flow  of  computations  in  adaptivity  with  a  go-back  procedure. 
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In  the  go-back  procedure,  two  different  meshes  are  used  for  the  computations  of  any 
time  interval.  Thus  the  adaptive  computation  will  be  at  least  twice  as  expensive  as  a 
coarse  mesh  computation. 

The  refinement  procedure  involves  four  steps: 

1.  check  the  error  or  refinement  criteria  for  all  elements; 

2.  generate  new  nodes  and  elements  where  refinement  is  needed; 

3.  generate  new  starting  data  for  new  elements  and  nodes; 

4.  generate  new  data  for  boundary  conditions,  loads,  etc. 

In  first  step,  the  error  or  refinement  indicators  are  computed  by  the  external  program, 
according  to  the  equations  in  Section  2-4.  These  indicators  are  checked  against  a  user- 
input  tolerance,  and  any  element  which  exceeds  the  tolerance  is  listed  for  subdivision. 
Only  elements  which  have  not  exceeded  a  specified  number  of  levels  of  refinement  are 
refined.  In  the  refinement  procedure,  new  elements  and  nodes  are  generated.  The  new 
nodes  and  elements  are  numbered  according  to  the  following:  the  lower-left  element 
retains  the  same  number  as  the  original  element,  and  the  other  three  new  elements  are 
added  to  the  end  of  the  element  list;  the  new  nodes  are  always  added  to  the  end  of  the 
node  list.  Nodes  are  classified  as  free  nodes  or  slave  nodes  depending  on  the 
refinement  of  the  contiguous  elements. 

Slave  nodes  which  are  generated  by  refinement  are  treated  similarly  to  tielines  in 
DYNA2D.  The  actual  DYNA2D  routines  could  not  be  used  for  the  slave  nodes  so  a 
modified  routine  was  written  which  implements  Equations  (2.10)  and  (2.14-2.15). 

The  starting  values  of  stresses  and  other  state  variables  must  be  set  in  the  elements 
which  have  been  generated  by  the  refinement,  and  the  velocities  and  displacements  are 
initiated  for  the  new  nodes.  The  new  starting  data  are  obtained  by  the  REZONE 
algorithm  in  DYNA2D,  which  sets  new  values  of  element  and  nodal  variables  according 
to  a  least  square  projection.  The  descendant  elements  are  assigned  the  same  material 
properties  as  the  parent  element. 
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The  fourth  step  requires  the  following  data  to  be  provided  so  that  it  reflects  the  addition 
of  the  new  nodes  to  the  model: 

1.  boundary  conditions 

2.  nodal  loads,  i.e.  f 

3.  slideline  data 

This  fourth  step  is  quite  difficult,  and  it  has  led  to  problems  in  this  program  and  many 
other  adaptive  programs,  for  it  requires  an  interpretation  of  the  problem  statement  from 
the  discrete  data.  In  other  words,  it  is  necessary  to  interpret  the  data  for  the  initial  mesh 
in  order  to  develop  the  data  for  the  new  nodes. 

For  the  boimdary  conditions,  the  following  rule  is  used.  The  original  nodes  are  called 
the  edge  parent  nodes.  The  botmdary  condition  of  the  new  node  is  set  according  to  the 
boundary  conditions  of  the  two  parent  nodes,  if  the  parent  nodes  have  the  same 
boundary  condition.  If  the  two  parent  nodes  have  different  boimdary  conditions,  but 
the  constraints  on  one  node  are  less  restrictive,  then  the  less  restrictive  boundary 
conditions  are  assigned  to  the  new  node.  For  example,  if  the  boundary  conditions  for 
parent  nodes  A  and  B  are: 

node  A:  x-component  fixed,  y-component  fixed; 
node  B:  x-component  fixed; 

then  the  x-component  is  fixed  for  the  new  node  between  nodes  A  and  B.  If  the 
boundary  conditions  of  the  nodes  are  such  that  the  conditions  on  one  parent  node  are 
not  a  subset  of  the  conditions  of  the  other  parent  node,  an  error  is  flagged. 

For  external  loads  specified  by  a  pressure  or  shear  load  card  (page  157  in  the  DYNA2D 
manual)  the  new  nodes  are  added  to  the  data.  Therefore,  the  pressure  and  shear  loads 
are  generated  exactly  as  in  a  regular  rim. 

A  similar  procedure  was  tried  for  the  slideline  data.  However,  we  have  not  been  able  to 
successfully  continue  a  run  with  slidelines  in  which  adaptivity  introduced  new  nodes 
along  the  slidelines.  This  option  requires  further  work. 

The  program  for  these  functions  is  in  a  block  of  subroutines  with  a  root  LP.  These 
subroutines  are  called  from  DYNA2D.  Two  ways  of  making  adaptive  nms  have  been 
included: 
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1.  automatic  adaptivity,  which  completes  a  run  using  the  go-back  procedure 
previously; 

2.  interactive  adaptivity  procediu'e,  in  which  the  user  can  initiate  adaptivity  at  any  time 
on  the  basis  on  run-time  plots  by  interrupting  the  run. 
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SECTIONS 

EXAMPLES  OF  H-ADAPTIVE  CALCULATIONS 
3.1  TAYLOR  IMPACT  PROBLEM. 

To  illustrate  the  performance  of  the  Ii-adaptive  procedure,  the  Taylor  rod  impact 
problem,  which  is  a  standard  DYNA2D  problem  was  solved.  The  problem  was  solved 
in  three  ways:  using  a  coarse  mesh  of  3x10  elements,  using  a  fine  mesh  of  12x10 
elements,  and  using  the  li-adaptive  procedure  starting  from  the  coarse  mesh.  The 
problem  is  axisymmetric,  so  axisymmetric  elements  were  used.  In  the  /i-adaptive 
simulation,  the  MLS  error  criterion  was  used  with  a  refinement  threshold  of  0.002  and 
two  levels  of  adaptivity  were  permitted.  The  total  simulation  time  was  80ms.  Error 
checks  and  adaptivity  were  performed  every  Sms,  so  a  changes  were  made  in  the  mesh 
16  times. 

The  rod  is  a  copper  rod  with  a  Yoxmg's  modulus  117GPa  and  a  Poisson's  ratio  of  0.4. 
The  material  is  isotropic-elastic-plastic  with  a  yield  stress  of  400  MPa  and  a  target 
modulus  of  100  MPa.  The  initial  velocity  of  the  rod  was  227m/s,  and  the  bottom  nodes 
were  constrained  in  the  z-direction  at  the  start  of  the  simulation  to  model  an  impact 
with  a  rigid  wall. 

The  evolution  of  the  adaptive  mesh  is  shown  in  Figures  3-1  and  3-2.  It  can  be  seen  that 
the  refinement  slowly  progresses  upward  from  the  impacted  surface,  but  that  the  top  of 
the  rod  is  never  refined.  For  comparison,  the  deformed  mesh  for  the  fine-mesh  solution 
is  shown  in  Figure  3-3.  Comparison  of  the  two  meshes  reveals  good  agreement. 

The  displacement  at  the  centerpoint  of  the  free  end  and  midpoint  of  the  rod  are 
compared  for  the  three  different  meshes  in  Figure  3-4.  It  can  be  seen  that  the  /z-adaptive 
solution  agrees  exactly  with  the  fine  mesh  solution,  whereas  the  coarse  mesh  solution 
differs  moderately  from  the  two.  The  velocities  for  these  two  points  are  compared  in 
Figure  3-5. 

Again,  the  adaptive  results  compare  very  well  with  the  fine  mesh  results.  At  several 
points  in  time,  the  //-adaptive  solution  exhibits  moderate  oscillations.  In  some  instances, 
these  appear  to  be  associated  with  the  refinement  of  the  elements  near  the  node  where 
the  velocity  is  output. 
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Figure  3-1.  Deformed  meshes  for  the  Taylor  impact  problem  for  /i-adaptive  solution  at 
5ps  intervals  from  5jis  to  40|js. 


’mesh’ 


in 

CO 


19 


intervals  from  45^5  to  80ps. 
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Figure  3-3.  Deformed  meshes 
to  80|is. 


J 

Figure  3-5.  Axial  velocities  of  free-end  point  and  midpoint  on  the  axis  of  symmetry 
for  Taylor  rod  problem;  time  is  in  (is. 
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Figxire  3-6  shows  the  axial  stress  a„  at  the  midpoint  of  the  rod  on  the  axis  of  S5nnmetry. 
The  /i-adaptive  agrees  very  well  with  the  fine  mesh  solution,  particularly  at  later  times, 
but  file  li-adaptive  solution  exhibits  severe  noise  at  about  40ms.  This  corresponds  to  the 
time  when  refinement  occurs  in  the  area  of  the  node  for  which  the  velocity  is  given.  The 
noise  decays  quite  rapidly  in  time  and  appears  to  have  no  permanent  effect  on  the 
solution. 

The  fine  mesh  solution  requires  55  minutes  of  computer  time,  tiie  /i-adaptive  solution 
required  35  minutes.  A  SUN-SPARC 1  work  station  was  used.  The  savings  in  computer 
time  are  moderate  but  nevertheless  significant.  It  should  be  noted  fiiat  substantially 
more  time  could  be  saved  if  the  DYNA2D  program  had  a  subcycling  capability  so  that 
different  time  steps  could  be  used  in  different  parts  of  the  mesh.  DYNA2D  uses  a  single 
time  step  for  the  entire  mesh,  so  that  when  the  fi-adaptivity  introduces  smaller  elements, 
the  integration  time  step  for  the  entire  mesh  is  needed. 
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Figure  3.6  The  axial  stress  Ozz  at  the  midpoint  of  the  axis  of  the  rod;  time  is  in  ps. 
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3.2  BEAM  PROBLEMS. 


The  response  of  a  cantilever  beam  was  simulated  for  an  elastic  material  and  an  elastic- 
plastic  material.  The  problem  definition  is  given  in  Figure  3-7,  all  units  are  non- 
dimensional.  A  vertical  load  of  magnitude  40  which  remains  in  the  vertical  direction  is 
applied  as  a  step  function  in  time  at  the  free  end;  the  load  is  a  traction  with  a  parabolic 
distribution  which  is  maximum  at  the  centerline  and  vanishes  at  the  top  and  bottom. 
The  material  parameters  are  as  follows:  Yotmg's  modulus=10,000,  Poisson's  ratio=300, 
yield  stress=300,  the  plastic  tangent  modulus=100,  and  the  density  =1;  all  units  are  non- 
dimensional. 

The  problem  was  solved  with  three  meshes:  a  coarse  mesh  with  2  elements  in  the 
vertical  direction,  10  elements  in  the  horizontal  direction;  a  fine  mesh  with  8  elements  in 
the  vertical  direction  and  40  elements  in  the  horizontal  direction;  an  Iz-adaptive  solution 
starting  with  the  coarse  mesh  and  two  levels  of  refinement.  For  the  adaptive  solution, 
the  projected  Green  strain  error  indicator  was  used  with  a  threshold  of  0.005  for 
refinement. 

The  time-histories  of  the  deflection  of  the  end-point  for  the  three  runs  for  the  elastic 
beam  is  shown  in  Figure  3-8;  for  the  elastic  solution,  no  yielding  was  allowed.  It  can  be 
seen  that  the  adaptive  solution  agrees  reasonably  well  with  the  fine-mesh  solution, 
although  the  agreement  between  the  /i-adaptive  and  fine  mesh  solutions  is  not  perfect. 
Some  of  the  discrepancy  may  be  caused  by  the  differences  between  the  loading  in  the 
two  solutions,  which  results  from  the  fact  that  in  the  fine  mesh  run,  the  load  was 
distributed  over  five  nodes  at  the  end,  whereas  for  the  coarse  mesh  solution  the  load  is 
distributed  over  3  nodes. 

The  time  histories  of  the  end-point  deflections  for  the  three  runs  with  the  elastic-plastic 
beam  are  shown  in  Figure  3-9.  Again,  there  are  moderate  differences  between  the  h- 
adaptive  solution  and  the  fine  mesh  solution,  but  the  ft-adaptive  solution  agrees  much 
better  with  the  fine  mesh  solution.  The  deformed  meshes  are  shown  in  Figure  3-10.  As 
can  be  seen,  the  refinement  is  concentrated  at  the  supported  edge  where  the  strains  and 
strain  gradients  are  maximum. 
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Figure  3-8.  Displacement  of  free  end  of  elastic  beam  as  a  function  of  time. 
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Figure  3-10.  Deformed  meshes  for  the  elastic-plastic  beam  at  four  stages  of 
deformation. 


SECTION  4 
CONCLUSIONS 


An  adaptive  procedure  has  been  implemented  in  the  DYNA2D  program  for  the 
nonlinear  transient  analysis  of  solids.  A  procedure  was  developed  so  that  the  program 
for  implementing  ft-adaptivity  is  totally  separate  from  the  DYNA2D  program.  This 
enables  the  procedure  to  be  used  with  other  versions  of  the  DYNA2D  program,  such  as 
subsequent  releases  or  modified  versions  developed  in-house  by  various  organizations. 
The  new  procedure  utilizes  the  rezone  and  remap  algorithms  in  DYNA2D  in 
conjimction  with  the  restart  options.  Thus  the  programming  for  adaptivity  deals 
primarily  with  the  implementation  of  the  error  and  refinement  criteria  and  with  the 
selection  of  elements  to  be  subdivided  and  is  completely  separate  from  DYNA2D. 

In  the  ^-adaptive  procedxire  developed  here,  the  mesh  can  be  refined  to  several  levels  by 
subdividing  quadrilaterals  into  four  quadrilaterals.  The  error  and  refinement  criteria 
select  the  elements  which  are  to  be  refined.  Several  error  and  refinement  indicators 
were  incorporated  in  the  program:  a  projection-based  error  indicator  using  the  Green 
strain  and  refinement  indicators  based  on  effective  stress  and  effective  plastic  strain.  In 
preliminary  studies,  we  also  tested  error  indicators  based  on  projections  on  the  stress. 
We  found  that  error  indicators  were  not  as  useful  for  the  class  of  nonlinear  materials 
which  were  studied,  elastic-plastic  type  materials,  because  the  stress  in  the  plastic  zone 
is  spatially  quite  uniform  even  when  there  are  large  variations  in  the  strains.  In  many 
situations,  a  criterion  for  refinement  based  on  the  magnitude  of  the  effective  plastic 
strain  works  quite  well,  because  refinement  is  usually  needed  in  the  areas  of  most 
severe  nonlinearities  and  these  are  the  areas  where  the  effective  plastic  strain  is  very 
large.  The  projection  criterion  based  on  Green  strain  was  quite  effective  in  refining  in 
areas  of  plastic  response. 

Comparisons  were  made  between  coarse  mesh,  fine  mesh  and  /z-adaptive  solutions  for 
three  problems:  the  Taylor  impact  problem,  an  elastic  beam  problem  and  an  elastic- 
plastic  beam  problem.  The  results  show  that  the  ft-adaptive  solutions  compare  quite 
well  with  the  fine  mesh  solution  even  though  the  li-adaptive  solution  requires 
substantially  fewer  elements.  It  was  found  that  the  moving  least  square  projection 
indicator  based  on  Green  strain  developed  here  is  quite  effective  in  selecting  the 
elements  that  need  to  be  refined.  The  running  time  for  the  /z-adaptive  solutions  is 
approximately  30-40%  less  than  for  the  fine  mesh  solutions.  This  is  a  substantial  saving 
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but  has  not  been  achieved,  for  example  in  li-adaptive  procedures  which  are  directly 
incorporated  in  a  program.  For  comparable  problems,  fine  mesh  solutions  generally 
took  two  and  a  half  times  as  long  as  an  adaptive  solution. 

There  are  two  reasons  for  the  failure  to  achieve  larger  savings  of  computation  time: 

1.  the  use  of  an  external  program  to  implement  adaptivity  involves  the  writing  of 
restart  files  and  the  transfer  between  the  two  programs  which  is  not  necessary  when  h- 
adaptivity  is  built  into  a  program; 

2.  the  DYNA2D  program  does  not  have  a  subcycling,  in  which  different  time  steps  can 
be  used  in  different  parts  of  the  mesh,  so  that  the  smaller  elements  can  be  integrated  by 
a  smaller  time  step,  while  the  remainder  of  the  mesh  is  integrated  with  larger  time  step. 

The  absence  of  subcycling  is  particularly  costly  m  li-adaptive  procedures  because  every 
level  of  refinement  decreases  the  stable  time  step  by  a  factor  of  two.  Therefore  it  is 
important  that  explicit  structural  dynamics  program  with  /i-adaptivity  have  a 
subcycling  feature.  It  is  recommended  that  subcycling  be  incorporated  in  the  DYNA2D 
codes  for  it  not  only  improves  the  performance  of  the  problems  in  /i-adaptive  problems 
but  can  also  lead  to  substantial  computational  savings  in  standard  solution  procedures. 

Another  difficulty  which  we  encoimtered  was  in  defining  boundary  conditions,  loads 
and  sliding  interface  data  for  meshes  subsequent  to  li-adaptivity.  With  the  structure  of 
current  programs  such  as  DYNA2D,  the  problem  definition  cind  consequently  the 
boimdary  conditions,  loads,  and  other  data  are  defined  discretely  and  no  higher  level 
description  of  the  problem  is  available  which  could  provide  the  primitive  data  is 
available.  This  is  a  particular  difficulty  in  prescribing  botmdary  conditions  to 
descendant  nodes  which  are  generated  in  the  /i-adaptive  procedure,  because  the 
boimdary  conditions  on  the  new  nodes  must  be  developed  from  simple  rules.  Although 
in  most  cases,  simple  rules  such  as  the  ones  devised  here  are  able  to  assign  the  right 
botmdary  conditions,  the  procedure  has  a  degree  of  uncertainty  which  is  not  desirable 
in  production  codes.  The  usability  of  programs  such  as  DYNA2D  could  be  enhanced  if 
the  descriptive  data  could  be  input  in  a  much  more  general  fashion  for  it  would 
eliminate  a  substantial  burden  for  the  users. 

The  procedure  developed  here  is  a  useful  model  for  introducing  the  li-adaptivity  into 
other  general  purpose  programs.  The  addition  of  adaptivity  by  directly  modifying  a 
program  is  usually  extremely  difficult  because  of  the  large  varieties  of  functionalities 
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which  characterize  these  programs.  When  internal  extensive  modifications  are  made,  it 
is  difficult  not  to  disrupt  these  functionalities.  On  the  other  hand,  restart  and  rezoning 
should  be  a  feature  of  any  useful  nonlinear  solid  mechanics  program.  By  utilizing  these 
features  as  described  here,  ^-adaptivity  can  be  implemented  in  a  reasonably 
straightforward  manner,  altinough  the  efficiency  will  not  match  that  of  a  program 
developed  completely  with  ft-adaptivity. 
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